home *** CD-ROM | disk | FTP | other *** search
/ Delphi Magazine Collection 2001 / Delphi Magazine Collection 20001 (2001).iso / DISKS / Issue69 / OpenGL / Geometry.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  2000-12-30  |  71.9 KB  |  2,073 lines
  1. // 22/03/00 - Egg - Added MakeShadowMatrix (adapted from "OpenGL SuperBible" book)
  2. // 21/03/00 - Egg - Removed PWordArray (was a SysUtils's duplicate)
  3. // 06/02/00 - Egg - Added VectorEquals
  4. // 05/02/00 - Egg - Added some "const", more still needed
  5. //                  Added overloads for some of the MakeXXXVector funcs
  6. //                  Added homogeneous vector consts, VectorSpacing
  7. // I've added some perf enhancements, but basicly, dont use these funcs if speed
  8. // is an issue (especially the matrix stuff), they are still too memory intensive.
  9. // Prefer direct vector operations.
  10. //
  11. unit Geometry;
  12.  
  13. // This unit contains many needed types, functions and procedures for
  14. // quaternion, vector and matrix arithmetics. It is specifically designed
  15. // for geometric calculations within R3 (affine vector space)
  16. // and R4 (homogeneous vector space).
  17. //
  18. // Note: The terms 'affine' or 'affine coordinates' are not really correct here
  19. //       because an 'affine transformation' describes generally a transformation which leads
  20. //       to a uniquely solvable system of equations and has nothing to do with the dimensionality
  21. //       of a vector. One could use 'projective coordinates' but this is also not really correct
  22. //       and since I haven't found a better name (or even any correct one), 'affine' is as good
  23. //       as any other one.
  24. //
  25. // Identifiers containing no dimensionality (like affine or homogeneous)
  26. // and no datatype (integer..extended) are supposed as R4 representation
  27. // with 'single' floating point type (examples are TVector, TMatrix,
  28. // and TQuaternion). The default data type is 'single' ('GLFloat' for OpenGL)
  29. // and used in all routines (except conversions and trigonometric functions).
  30. //
  31. // Routines with an open array as argument can either take Func([1,2,3,4,..]) or Func(Vect).
  32. // The latter is prefered, since no extra stack operations is required.
  33. // Note: Be careful while passing open array elements! If you pass more elements
  34. //       than there's room in the result the behaviour will be unpredictable.
  35. //
  36. // If not otherwise stated, all angles are given in radians
  37. // (instead of degrees). Use RadToDeg or DegToRad to convert between them.
  38. //
  39. // Geometry.pas was assembled from different sources (like GraphicGems)
  40. // and relevant books or based on self written code, respectivly.
  41. //
  42. // Note: Some aspects need to be considered when using Delphi and pure
  43. //       assembler code. Delphi esnures that the direction flag is always
  44. //       cleared while entering a function and expects it cleared on return.
  45. //       This is in particular important in routines with (CPU) string commands (MOVSD etc.)
  46. //       The registers EDI, ESI and EBX (as well as the stack management
  47. //       registers EBP and ESP) must not be changed! EAX, ECX and EDX are
  48. //       freely available and mostly used for parameter.
  49. //
  50. // Version 2.5
  51. // last change : 04. January 2000
  52. //
  53. // (c) Copyright 1999, Dipl. Ing. Mike Lischke (public@lischke-online.de)
  54.  
  55. interface
  56.  
  57. type
  58.   // data types needed for 3D graphics calculation,
  59.   // included are 'C like' aliases for each type (to be
  60.   // conformal with OpenGL types)
  61.  
  62.   PByte = ^Byte;
  63.   PWord = ^Word;
  64.   PInteger = ^Integer;
  65.   PFloat = ^Single;
  66.   PDouble = ^Double;
  67.   PExtended = ^Extended;
  68.   PPointer = ^Pointer;
  69.  
  70.   // types to specify continous streams of a specific type
  71.   // switch off range checking to access values beyond the limits 
  72.   PByteVector = ^TByteVector;
  73.   PByteArray = PByteVector;
  74.   TByteVector = array[0..0] of Byte;
  75.  
  76.   PWordVector = ^TWordVector;
  77.   TWordVector = array[0..0] of Word;
  78.  
  79.   PIntegerVector = ^TIntegerVector;
  80.   PIntegerArray = PIntegerVector;
  81.   TIntegerVector = array[0..0] of Integer;
  82.  
  83.   PFloatVector = ^TFloatVector;
  84.   PFloatArray = PFloatVector;
  85.   TFloatVector = array[0..0] of Single;
  86.  
  87.   PDoubleVector = ^TDoubleVector;
  88.   PDoubleArray = PDoubleVector;
  89.   TDoubleVector = array[0..0] of Double;
  90.  
  91.   PExtendedVector = ^TExtendedVector;
  92.   PExtendedArray = PExtendedVector;
  93.   TExtendedVector = array[0..0] of Extended;
  94.  
  95.   PPointerVector = ^TPointerVector;
  96.   PPointerArray = PPointerVector;
  97.   TPointerVector = array[0..0] of Pointer;
  98.  
  99.   PCardinalVector = ^TCardinalVector;
  100.   PCardinalArray = PCardinalVector;
  101.   TCardinalVector = array[0..0] of Cardinal;
  102.  
  103.   // common vector and matrix types with predefined limits
  104.   // indices correspond like: x -> 0
  105.   //                          y -> 1
  106.   //                          z -> 2
  107.   //                          w -> 3
  108.  
  109.   PHomogeneousByteVector = ^THomogeneousByteVector;
  110.   THomogeneousByteVector = array[0..3] of Byte;
  111.   TVector4b = THomogeneousByteVector;
  112.  
  113.   PHomogeneousWordVector = ^THomogeneousWordVector;
  114.   THomogeneousWordVector = array[0..3] of Word;
  115.   TVector4w = THomogeneousWordVector;
  116.  
  117.   PHomogeneousIntVector = ^THomogeneousIntVector;
  118.   THomogeneousIntVector = array[0..3] of Integer;
  119.   TVector4i = THomogeneousIntVector;
  120.  
  121.   PHomogeneousFltVector = ^THomogeneousFltVector;
  122.   THomogeneousFltVector = array[0..3] of Single;
  123.   TVector4f = THomogeneousFltVector;
  124.  
  125.   PHomogeneousDblVector = ^THomogeneousDblVector;
  126.   THomogeneousDblVector = array[0..3] of Double;
  127.   TVector4d = THomogeneousDblVector;
  128.  
  129.   PHomogeneousExtVector = ^THomogeneousExtVector;
  130.   THomogeneousExtVector = array[0..3] of Extended;
  131.   TVector4e = THomogeneousExtVector;
  132.  
  133.   PHomogeneousPtrVector = ^THomogeneousPtrVector;
  134.   THomogeneousPtrVector = array[0..3] of Pointer;
  135.   TVector4p = THomogeneousPtrVector;
  136.  
  137.   PAffineByteVector = ^TAffineByteVector;
  138.   TAffineByteVector = array[0..2] of Byte;
  139.   TVector3b = TAffineByteVector;
  140.  
  141.   PAffineWordVector = ^TAffineWordVector;
  142.   TAffineWordVector = array[0..2] of Word;
  143.   TVector3w = TAffineWordVector;
  144.  
  145.   PAffineIntVector = ^TAffineIntVector;
  146.   TAffineIntVector = array[0..2] of Integer;
  147.   TVector3i = TAffineIntVector;
  148.  
  149.   PAffineFltVector = ^TAffineFltVector;
  150.   TAffineFltVector = array[0..2] of Single;
  151.   TVector3f = TAffineFltVector;
  152.  
  153.   PAffineDblVector = ^TAffineDblVector;
  154.   TAffineDblVector = array[0..2] of Double;
  155.   TVector3d = TAffineDblVector;
  156.  
  157.   PAffineExtVector = ^TAffineExtVector;
  158.   TAffineExtVector = array[0..2] of Extended;
  159.   TVector3e = TAffineExtVector;
  160.  
  161.   PAffinePtrVector = ^TAffinePtrVector;
  162.   TAffinePtrVector = array[0..2] of Pointer;
  163.   TVector3p = TAffinePtrVector;
  164.  
  165.   // some simplified names
  166.   PVector = ^TVector;
  167.   TVector = THomogeneousFltVector;
  168.  
  169.   PHomogeneousVector = ^THomogeneousVector;
  170.   THomogeneousVector = THomogeneousFltVector;
  171.  
  172.   PAffineVector = ^TAffineVector;
  173.   TAffineVector = TAffineFltVector;
  174.  
  175.   // arrays of vectors
  176.   PVectorArray = ^TVectorArray;
  177.   TVectorArray = array[0..0] of TAffineVector;
  178.  
  179.   // matrices
  180.   THomogeneousByteMatrix = array[0..3] of THomogeneousByteVector;
  181.   TMatrix4b = THomogeneousByteMatrix;
  182.  
  183.   THomogeneousWordMatrix = array[0..3] of THomogeneousWordVector;
  184.   TMatrix4w = THomogeneousWordMatrix;
  185.  
  186.   THomogeneousIntMatrix = array[0..3] of THomogeneousIntVector;
  187.   TMatrix4i = THomogeneousIntMatrix;
  188.  
  189.   THomogeneousFltMatrix  = array[0..3] of THomogeneousFltVector;
  190.   TMatrix4f = THomogeneousFltMatrix;
  191.  
  192.   THomogeneousDblMatrix = array[0..3] of THomogeneousDblVector;
  193.   TMatrix4d = THomogeneousDblMatrix;
  194.  
  195.   THomogeneousExtMatrix = array[0..3] of THomogeneousExtVector;
  196.   TMatrix4e = THomogeneousExtMatrix;
  197.  
  198.   TAffineByteMatrix = array[0..2] of TAffineByteVector;
  199.   TMatrix3b = TAffineByteMatrix;
  200.  
  201.   TAffineWordMatrix = array[0..2] of TAffineWordVector;
  202.   TMatrix3w = TAffineWordMatrix;
  203.  
  204.   TAffineIntMatrix = array[0..2] of TAffineIntVector;
  205.   TMatrix3i = TAffineIntMatrix;
  206.  
  207.   TAffineFltMatrix = array[0..2] of TAffineFltVector;
  208.   TMatrix3f = TAffineFltMatrix;
  209.  
  210.   TAffineDblMatrix = array[0..2] of TAffineDblVector;
  211.   TMatrix3d = TAffineDblMatrix;
  212.  
  213.   TAffineExtMatrix = array[0..2] of TAffineExtVector;
  214.   TMatrix3e = TAffineExtMatrix;
  215.  
  216.   // some simplified names
  217.   PMatrix = ^TMatrix;
  218.   TMatrix = THomogeneousFltMatrix;
  219.  
  220.   PHomogeneousMatrix = ^THomogeneousMatrix;
  221.   THomogeneousMatrix = THomogeneousFltMatrix;
  222.  
  223.   PAffineMatrix = ^TAffineMatrix;
  224.   TAffineMatrix = TAffineFltMatrix;
  225.  
  226.   // q = ([x, y, z], w)
  227.   TQuaternion = record
  228.     case Integer of
  229.       0:
  230.         (ImagPart: TAffineVector;
  231.          RealPart: Single);
  232.       1:
  233.         (Vector: TVector4f);
  234.   end;
  235.  
  236.   TRectangle = record
  237.     Left,
  238.     Top,
  239.     Width,
  240.     Height: Integer;
  241.   end;
  242.  
  243.   TTransType = (ttScaleX, ttScaleY, ttScaleZ,
  244.                 ttShearXY, ttShearXZ, ttShearYZ,
  245.                 ttRotateX, ttRotateY, ttRotateZ,
  246.                 ttTranslateX, ttTranslateY, ttTranslateZ,
  247.                 ttPerspectiveX, ttPerspectiveY, ttPerspectiveZ, ttPerspectiveW);
  248.  
  249.   // used to describe a sequence of transformations in following order:
  250.   // [Sx][Sy][Sz][ShearXY][ShearXZ][ShearZY][Rx][Ry][Rz][Tx][Ty][Tz][P(x,y,z,w)]
  251.   // constants are declared for easier access (see MatrixDecompose below)
  252.   TTransformations  = array [TTransType] of Single;
  253.  
  254.     
  255. const
  256.   // useful constants
  257.  
  258.   // standard vectors
  259.   XVector :    TAffineVector = (1, 0, 0);
  260.   YVector :    TAffineVector = (0, 1, 0);
  261.   ZVector :    TAffineVector = (0, 0, 1);
  262.   XYZVector :  TAffineVector = (1, 1, 1);
  263.   NullVector : TAffineVector = (0, 0, 0);
  264.   // standard homogeneous vectors
  265.   XHmgVector : THomogeneousVector = (1, 0, 0, 0);
  266.   YHmgVector : THomogeneousVector = (0, 1, 0, 0);
  267.   ZHmgVector : THomogeneousVector = (0, 0, 1, 0);
  268.   WHmgVector : THomogeneousVector = (0, 0, 0, 1);
  269.   XYZHmgVector  : THomogeneousVector = (1, 1, 1, 0);
  270.   XYZWHmgVector : THomogeneousVector = (1, 1, 1, 1);
  271.   NullHmgVector : THomogeneousVector = (0, 0, 0, 0);
  272.   // standard homogeneous points
  273.   XHmgPoint :  THomogeneousVector = (1, 0, 0, 1);
  274.   YHmgPoint :  THomogeneousVector = (0, 1, 0, 1);
  275.   ZHmgPoint :  THomogeneousVector = (0, 0, 1, 1);
  276.   WHmgPoint :  THomogeneousVector = (0, 0, 0, 1);
  277.   NullHmgPoint : THomogeneousVector = (0, 0, 0, 1);
  278.  
  279.   IdentityMatrix: TMatrix = ((1, 0, 0, 0),
  280.                              (0, 1, 0, 0),
  281.                              (0, 0, 1, 0),
  282.                              (0, 0, 0, 1));
  283.   EmptyMatrix: TMatrix = ((0, 0, 0, 0),
  284.                           (0, 0, 0, 0),
  285.                           (0, 0, 0, 0),
  286.                           (0, 0, 0, 0));
  287.   // some very small numbers
  288.   EPSILON  = 1e-100;
  289.   EPSILON2 = 1e-50;
  290.  
  291. //----------------------------------------------------------------------------------------------------------------------
  292.  
  293. // vector functions
  294. function  VectorAdd(const V1, V2: TVector): TVector;
  295. //: Returns the sum of two vectors
  296. function  VectorAffineAdd(const V1, V2: TAffineVector): TAffineVector;
  297. //: Makes a linear combination of two vectors and return the result
  298. function  VectorAffineCombine(const V1, V2: TAffineVector; F1, F2: Single): TAffineVector;
  299. //: Makes a linear combination of three vectors and return the result
  300. function  VectorAffineCombine3(const V1, V2, V3: TAffineVector; F1, F2, F3: Single): TAffineVector;
  301. function  VectorAffineDotProduct(const V1, V2: TAffineVector): Single;
  302. function  VectorAffineLerp(const V1, V2: TAffineVector; t: Single): TAffineVector;
  303. function  VectorAffineSubtract(const V1, V2: TAffineVector): TAffineVector;
  304. function  VectorAngle(const V1, V2: TAffineVector): Single;
  305. function  VectorCombine(const V1, V2: TVector; F1, F2: Single): TVector;
  306. //: Calculates the cross product between vector 1 and 2
  307. function  VectorCrossProduct(const V1, V2: TAffineVector): TAffineVector; overload;
  308. function  VectorCrossProduct(const V1, V2: TVector): TVector; overload;
  309. function  VectorDotProduct(const V1, V2: TVector): Single;
  310. function  VectorLength(V: array of Single): Single;
  311. function  VectorLerp(const V1, V2: TVector; t: Single): TVector;
  312. procedure VectorNegate(V: array of Single);
  313. //: Calculates norm of a vector which is defined as norm = x * x + y * y + ...
  314. function  VectorNorm(V: array of Single): Single;
  315. //: Transforms a vector to unit length and return length
  316. function  VectorNormalize(V: array of Single): Single;
  317. {: Calculates a vector perpendicular to N.<p>
  318.    N is assumed to be of unit length,  subtract out any component parallel to N }
  319. function  VectorPerpendicular(const V, N: TAffineVector): TAffineVector;
  320. function  VectorReflect(const V, N: TAffineVector): TAffineVector;
  321. //: Rotates Vector about Axis with Angle radiants
  322. procedure VectorRotate(var Vector: TVector4f; const Axis: TVector3f; Angle: Single);
  323. //: Returns a vector scaled by a factor
  324. procedure VectorScale(V: array of Single; Factor: Single);
  325. //: Returns V1-V2
  326. function  VectorSubtract(const V1, V2: TVector): TVector;
  327. {: Calculates Abs(v1[x]-v2[x])+Abs(v1[y]-v2[y])+ etc.<p> }
  328. function  VectorSpacing(const V1, V2: TVector): Single;
  329. //: True if all components are equal.
  330. function  VectorEquals(const V1, V2: TVector) : Boolean; overload;
  331. //: True if all components are equal.
  332. function  VectorEquals(const V1, V2: TAffineVector) : Boolean; overload;
  333.  
  334. // matrix functions
  335. function  CreateRotationMatrixX(Sine, Cosine: Single): TMatrix;
  336. function  CreateRotationMatrixY(Sine, Cosine: Single): TMatrix;
  337. function  CreateRotationMatrixZ(Sine, Cosine: Single): TMatrix;
  338. function  CreateScaleMatrix(const V: TAffineVector): TMatrix;
  339. function  CreateTranslationMatrix(const V: TVector): TMatrix;
  340. procedure MatrixAdjoint(var M: TMatrix);
  341. function  MatrixAffineDeterminant(const M: TAffineMatrix): Single;
  342. procedure MatrixAffineTranspose(var M: TAffineMatrix);
  343. function  MatrixDeterminant(const M: TMatrix): Single;
  344. procedure MatrixInvert(var M: TMatrix);
  345. function  MatrixMultiply(const M1, M2: TMatrix): TMatrix;
  346. procedure MatrixScale(var M: TMatrix; Factor: Single);
  347. procedure MatrixTranspose(var M: TMatrix);
  348.  
  349. // quaternion functions
  350. function  QuaternionConjugate(const Q: TQuaternion): TQuaternion;
  351. function  QuaternionFromPoints(const V1, V2: TAffineVector): TQuaternion;
  352. function  QuaternionMultiply(const qL, qR: TQuaternion): TQuaternion;
  353. function  QuaternionSlerp(const QStart, QEnd: TQuaternion; Spin: Integer; t: Single): TQuaternion;
  354. function  QuaternionToMatrix(const Q: TQuaternion): TMatrix;
  355. procedure QuaternionToPoints(const Q: TQuaternion; var ArcFrom, ArcTo: TAffineVector);
  356.  
  357. // mixed functions
  358. function  ConvertRotation(const Angles: TAffineVector): TVector;
  359. function  CreateRotationMatrix(const Axis: TVector3f; Angle: Single): TMatrix;
  360. function  MatrixDecompose(const M: TMatrix; var Tran: TTransformations): Boolean;
  361. //: Transforms an affine vector by multiplying it with a matrix
  362. function  VectorAffineTransform(const V: TAffineVector; const M: TAffineMatrix): TAffineVector;
  363. //: Transforms a homogeneous vector by multiplying it with a matrix
  364. function  VectorTransform(const V: TVector4f; const M: TMatrix): TVector4f; overload;
  365. //: Transforms an affine vector by multiplying it with a (homogeneous) matrix
  366. function  VectorTransform(const V: TVector3f; const M: TMatrix): TVector3f; overload;
  367.  
  368. // miscellaneous functions
  369. function  MakeAffineDblVector(V: array of Double): TAffineDblVector;
  370. function  MakeDblVector(V: array of Double): THomogeneousDblVector;
  371. function  MakeAffineVector(V: array of Single) : TAffineVector; overload
  372. function  MakeAffineVector(const x, y, z : Single) : TAffineVector; overload;
  373. function  MakeQuaternion(Imag: array of Single; Real: Single): TQuaternion;
  374. function  MakeVector(V: array of Single) : TVector; overload;
  375. function  MakeVector(const x, y, z: Single; w : Single = 0) : TVector; overload;
  376. function  MakePoint(const x, y, z: Single; w : Single = 1) : TVector; overload;
  377. function  PointInPolygon(xp, yp : array of Single; x, y: Single): Boolean;
  378. function  VectorAffineDblToFlt(const V: TAffineDblVector): TAffineVector;
  379. function  VectorDblToFlt(const V: THomogeneousDblVector): THomogeneousVector;
  380. function  VectorAffineFltToDbl(const V: TAffineVector): TAffineDblVector;
  381. function  VectorFltToDbl(const V: TVector): THomogeneousDblVector;
  382.  
  383. // trigonometric functions
  384. function  ArcCos(const X: Extended): Extended;
  385. function  ArcSin(const X: Extended): Extended;
  386. function  ArcTan2(const Y, X: Extended): Extended;
  387. function  CoTan(const X: Extended): Extended;
  388. function  DegToRad(const Degrees: Extended): Extended;
  389. function  RadToDeg(const Radians: Extended): Extended;
  390. //: Calculates sine and cosine from the given angle Theta
  391. procedure SinCos(const Theta: Extended; var Sin, Cos: Extended);
  392. function  Tan(const X: Extended): Extended;
  393.  
  394. // coordinate system manipulation functions
  395.  
  396. //: Rotates the given coordinate system (represented by the matrix) around its Y-axis
  397. function Turn(const Matrix: TMatrix; Angle: Single): TMatrix; overload;
  398. //: Rotates the given coordinate system (represented by the matrix) around MasterUp
  399. function Turn(const Matrix: TMatrix; const MasterUp: TAffineVector; Angle: Single): TMatrix; overload;
  400. //: Rotates the given coordinate system (represented by the matrix) around its X-axis
  401. function Pitch(const Matrix: TMatrix; Angle: Single): TMatrix; overload;
  402. //: Rotates the given coordinate system (represented by the matrix) around MasterRight
  403. function Pitch(const Matrix: TMatrix; const MasterRight: TAffineVector; Angle: Single): TMatrix; overload;
  404. //: Rotates the given coordinate system (represented by the matrix) around its Z-axis
  405. function Roll(const Matrix: TMatrix; Angle: Single): TMatrix; overload;
  406. //: Rotates the given coordinate system (represented by the matrix) around MasterDirection
  407. function Roll(const Matrix: TMatrix; const MasterDirection: TAffineVector; Angle: Single): TMatrix; overload;
  408.  
  409. // misc funcs
  410.  
  411. // Creates a shadow projection matrix out of the plane equation
  412. // coefficients and the position of the light. The return value is stored
  413. // in destMat[][]
  414. function MakeShadowMatrix(const planePoint, planeNormal, lightPos : TVector) : TMatrix;
  415.  
  416. //----------------------------------------------------------------------------------------------------------------------
  417.  
  418. implementation
  419.  
  420. const
  421.   // FPU status flags (high order byte)
  422.   C0 = 1;
  423.   C1 = 2;
  424.   C2 = 4;
  425.   C3 = $40;
  426.  
  427.   // to be used as descriptive indices
  428.   X = 0;
  429.   Y = 1;
  430.   Z = 2;
  431.   W = 3;
  432.  
  433. //----------------- trigonometric helper functions ---------------------------------------------------------------------
  434.  
  435. function DegToRad(const Degrees: Extended): Extended;
  436.  
  437. begin
  438.   Result := Degrees * (PI / 180);
  439. end;
  440.  
  441. //----------------------------------------------------------------------------------------------------------------------
  442.  
  443. function RadToDeg(const Radians: Extended): Extended;
  444.  
  445. begin
  446.   Result := Radians * (180 / PI);
  447. end;
  448.  
  449. // SinCos (Extended)
  450. //
  451. procedure SinCos(const Theta: Extended; var Sin, Cos: Extended); assembler; register;
  452. // EAX contains address of Sin
  453. // EDX contains address of Cos
  454. // Theta is passed over the stack
  455. asm
  456.    FLD  Theta
  457.    FSINCOS
  458.    FSTP TBYTE PTR [EDX]    // cosine
  459.    FSTP TBYTE PTR [EAX]    // sine
  460.    FWAIT
  461. end;
  462.  
  463. //----------------------------------------------------------------------------------------------------------------------
  464.  
  465. function ArcCos(const X: Extended): Extended;
  466.  
  467. begin
  468.   Result := ArcTan2(Sqrt(1 - X * X), X);
  469. end;
  470.  
  471. //----------------------------------------------------------------------------------------------------------------------
  472.  
  473. function ArcSin(const X: Extended): Extended;
  474.  
  475. begin
  476.   Result := ArcTan2(X, Sqrt(1 - X * X))
  477. end;
  478.  
  479. //----------------------------------------------------------------------------------------------------------------------
  480.  
  481. function ArcTan2(const Y, X: Extended): Extended;
  482.  
  483. asm
  484.               FLD  Y
  485.               FLD  X
  486.               FPATAN
  487.               FWAIT
  488. end;
  489.  
  490. //----------------------------------------------------------------------------------------------------------------------
  491.  
  492. function Tan(const X: Extended): Extended;
  493.  
  494. asm
  495.               FLD  X
  496.               FPTAN
  497.               FSTP ST(0)      // FPTAN pushes 1.0 after result
  498.               FWAIT
  499. end;
  500.  
  501. //----------------------------------------------------------------------------------------------------------------------
  502.  
  503. function CoTan(const X: Extended): Extended;
  504.  
  505. asm
  506.               FLD  X
  507.               FPTAN
  508.               FDIVRP
  509.               FWAIT
  510. end;
  511.  
  512. //----------------- miscellaneous vector functions ---------------------------------------------------------------------
  513.  
  514. function MakeAffineDblVector(V: array of Double): TAffineDblVector; assembler;
  515.  
  516. // creates a vector from given values
  517. // EAX contains address of V
  518. // ECX contains address to result vector
  519. // EDX contains highest index of V
  520.  
  521. asm
  522.               PUSH EDI
  523.               PUSH ESI
  524.               MOV EDI, ECX
  525.               MOV ESI, EAX
  526.               MOV ECX, EDX
  527.               ADD ECX, 2
  528.               REP MOVSD
  529.               POP ESI
  530.               POP EDI
  531. end;
  532.  
  533. //----------------------------------------------------------------------------------------------------------------------
  534.  
  535. function MakeDblVector(V: array of Double): THomogeneousDblVector; assembler;
  536.  
  537. // creates a vector from given values
  538. // EAX contains address of V
  539. // ECX contains address to result vector
  540. // EDX contains highest index of V
  541.  
  542. asm
  543.               PUSH EDI
  544.               PUSH ESI
  545.               MOV EDI, ECX
  546.               MOV ESI, EAX
  547.               MOV ECX, EDX
  548.               ADD ECX, 2
  549.               REP MOVSD
  550.               POP ESI
  551.               POP EDI
  552. end;
  553.  
  554. //----------------------------------------------------------------------------------------------------------------------
  555.  
  556. function MakeAffineVector(V: array of Single): TAffineVector; assembler;
  557. // creates a vector from given values
  558. // EAX contains address of V
  559. // ECX contains address to result vector
  560. // EDX contains highest index of V
  561. asm
  562.               PUSH EDI
  563.               PUSH ESI
  564.               MOV EDI, ECX
  565.               MOV ESI, EAX
  566.               MOV ECX, EDX
  567.               INC ECX
  568.               CMP ECX, 3
  569.               JB  @@1
  570.               MOV ECX, 3
  571. @@1:          REP MOVSD                     // copy given values
  572.               MOV ECX, 2
  573.               SUB ECX, EDX                   // determine missing entries
  574.               JS  @@Finish
  575.               XOR EAX, EAX
  576.               REP STOSD                     // set remaining fields to 0
  577. @@Finish:     POP ESI
  578.               POP EDI
  579. end;
  580.  
  581. // MakeAffineVector
  582. //
  583. function MakeAffineVector(const x, y, z : Single) : TAffineVector; overload;
  584. begin
  585.    Result[0]:=x;
  586.    Result[1]:=y;
  587.    Result[2]:=z;
  588. end;
  589.  
  590.  
  591. //----------------------------------------------------------------------------------------------------------------------
  592.  
  593. function MakeQuaternion(Imag: array of Single; Real: Single): TQuaternion; assembler;
  594.  
  595. // creates a quaternion from the given values
  596. // EAX contains address of Imag
  597. // ECX contains address to result vector
  598. // EDX contains highest index of Imag
  599. // Real part is passed on the stack
  600.  
  601. asm
  602.               PUSH EDI
  603.               PUSH ESI
  604.               MOV EDI, ECX
  605.               MOV ESI, EAX
  606.               MOV ECX, EDX
  607.               INC ECX
  608.               REP MOVSD
  609.               MOV EAX, [Real]
  610.               MOV [EDI], EAX
  611.               POP ESI
  612.               POP EDI
  613. end;
  614.  
  615. //----------------------------------------------------------------------------------------------------------------------
  616.  
  617. function MakeVector(V: array of Single): TVector; assembler;
  618. // creates a vector from given values
  619. // EAX contains address of V
  620. // ECX contains address to result vector
  621. // EDX contains highest index of V
  622. asm
  623.               PUSH EDI
  624.               PUSH ESI
  625.               MOV EDI, ECX
  626.               MOV ESI, EAX
  627.               MOV ECX, EDX
  628.               INC ECX
  629.               CMP ECX, 4
  630.               JB  @@1
  631.                   MOV ECX, 4
  632. @@1:          REP MOVSD                     // copy given values
  633.               MOV ECX, 3
  634.               SUB ECX, EDX                   // determine missing entries
  635.               JS  @@Finish
  636.               XOR EAX, EAX
  637.               REP STOSD                     // set remaining fields to 0
  638. @@Finish:     POP ESI
  639.               POP EDI
  640. end;
  641.  
  642. // MakeVector
  643. //
  644. function MakeVector(const x, y, z : Single; w : Single = 0) : TVector;
  645. begin
  646.     Result[0]:=x;
  647.     Result[1]:=y;
  648.     Result[2]:=z;
  649.     Result[3]:=w;
  650. end;
  651.  
  652. // MakePoint
  653. //
  654. function  MakePoint(const x, y, z: Single; w : Single = 1) : TVector;
  655. begin
  656.     Result[0]:=x;
  657.     Result[1]:=y;
  658.     Result[2]:=z;
  659.     Result[3]:=w;
  660. end;
  661.  
  662. //----------------------------------------------------------------------------------------------------------------------
  663.  
  664. function VectorLength(V: array of Single): Single; assembler;
  665.  
  666. // calculates the length of a vector following the equation: sqrt(x * x + y * y + ...)
  667. // Note: The parameter of this function is declared as open array. Thus
  668. // there's no restriction about the number of the components of the vector.
  669. //
  670. // EAX contains address of V
  671. // EDX contains the highest index of V
  672. // the result is returned in ST(0)
  673.  
  674. asm
  675.               FLDZ                           // initialize sum
  676. @@Loop:       FLD  DWORD PTR [EAX  +  4 * EDX] // load a component
  677.               FMUL ST, ST
  678.               FADDP
  679.               SUB  EDX, 1
  680.               JNL  @@Loop
  681.               FSQRT
  682. end;
  683.  
  684. //----------------------------------------------------------------------------------------------------------------------
  685.  
  686. function VectorAngle(const V1, V2: TAffineVector): Single; assembler;
  687.  
  688. // calculates the cosine of the angle between Vector1 and Vector2
  689. // Result = DotProduct(V1, V2) / (Length(V1) * Length(V2))
  690. //
  691. // EAX contains address of Vector1
  692. // EDX contains address of Vector2
  693.  
  694. asm
  695.               FLD DWORD PTR [EAX]           // V1[0]
  696.               FLD ST                        // double V1[0]
  697.               FMUL ST, ST                   // V1[0]^2 (prep. for divisor)
  698.               FLD DWORD PTR [EDX]           // V2[0]
  699.               FMUL ST(2), ST                // ST(2) := V1[0] * V2[0]
  700.               FMUL ST, ST                   // V2[0]^2 (prep. for divisor)
  701.               FLD DWORD PTR [EAX + 4]       // V1[1]
  702.               FLD ST                        // double V1[1]
  703.               FMUL ST, ST                   // ST(0) := V1[1]^2
  704.               FADDP ST(3), ST               // ST(2) := V1[0]^2 + V1[1] *  * 2
  705.               FLD DWORD PTR [EDX + 4]       // V2[1]
  706.               FMUL ST(1), ST                // ST(1) := V1[1] * V2[1]
  707.               FMUL ST, ST                   // ST(0) := V2[1]^2
  708.               FADDP ST(2), ST               // ST(1) := V2[0]^2 + V2[1]^2
  709.               FADDP ST(3), ST               // ST(2) := V1[0] * V2[0] + V1[1] * V2[1]
  710.               FLD DWORD PTR [EAX + 8]       // load V2[1]
  711.               FLD ST                        // same calcs go here
  712.               FMUL ST, ST                   // (compare above)
  713.               FADDP ST(3), ST
  714.               FLD DWORD PTR [EDX + 8]
  715.               FMUL ST(1), ST
  716.               FMUL ST, ST
  717.               FADDP ST(2), ST
  718.               FADDP ST(3), ST
  719.               FMULP                         // ST(0) := (V1[0]^2 + V1[1]^2 + V1[2]) *
  720.                                             //          (V2[0]^2 + V2[1]^2 + V2[2])
  721.               FSQRT                         // sqrt(ST(0))
  722.               FDIVP                         // ST(0) := Result := ST(1) / ST(0)
  723.   // the result is expected in ST(0), if it's invalid, an error is raised
  724. end;
  725.  
  726. // VectorNorm
  727. //
  728. function VectorNorm(V: array of Single): Single; assembler; register;
  729. // EAX contains address of V
  730. // EDX contains highest index in V
  731. // result is passed in ST(0)
  732. asm
  733.               FLDZ                           // initialize sum
  734. @@Loop:       FLD  DWORD PTR [EAX + 4 * EDX] // load a component
  735.               FMUL ST, ST                    // make square
  736.               FADDP                          // add previous calculated sum
  737.               SUB  EDX, 1
  738.               JNL  @@Loop
  739. end;
  740.  
  741. //----------------------------------------------------------------------------------------------------------------------
  742.  
  743. function VectorNormalize(V: array of Single): Single; assembler; register;
  744.  
  745. // EAX contains address of V
  746. // EDX contains the highest index in V
  747. // return former length of V in ST
  748.  
  749. asm
  750.               PUSH EBX
  751.               MOV ECX, EDX                  // save size of V
  752.               CALL VectorLength             // calculate length of vector
  753.               FTST                          // test if length = 0
  754.               MOV EBX, EAX                  // save parameter address
  755.               FSTSW AX                      // get test result
  756.               TEST AH, C3                   // check the test result
  757.               JNZ @@Finish
  758.               SUB EBX, 4                    // simplyfied address calculation
  759.               INC ECX
  760.               FLD1                          // calculate reciprocal of length
  761.               FDIV ST, ST(1)
  762. @@1:          FLD ST                        // double reciprocal
  763.               FMUL DWORD PTR [EBX + 4 * ECX] // scale component
  764.               WAIT
  765.               FSTP DWORD PTR [EBX + 4 * ECX] // store result
  766.               LOOP @@1
  767.               FSTP ST                       // remove reciprocal from FPU stack
  768. @@Finish:     POP EBX
  769. end;
  770.  
  771. //----------------------------------------------------------------------------------------------------------------------
  772.  
  773. function VectorAffineSubtract(const V1, V2: TAffineVector): TAffineVector; assembler; register;
  774.  
  775. // returns v1 minus v2
  776. // EAX contains address of V1
  777. // EDX contains address of V2
  778. // ECX contains address of the result
  779.  
  780. asm
  781.   {Result[X] := V1[X]-V2[X];
  782.   Result[Y] := V1[Y]-V2[Y];
  783.   Result[Z] := V1[Z]-V2[Z];}
  784.  
  785.                    FLD DWORD PTR [EAX]
  786.                    FSUB DWORD PTR [EDX]
  787.                    FSTP DWORD PTR [ECX]
  788.                    FLD DWORD PTR [EAX + 4]
  789.                    FSUB DWORD PTR [EDX + 4]
  790.                    FSTP DWORD PTR [ECX + 4]
  791.                    FLD DWORD PTR [EAX + 8]
  792.                    FSUB DWORD PTR [EDX + 8]
  793.                    FSTP DWORD PTR [ECX + 8]
  794. end;
  795.  
  796. //----------------------------------------------------------------------------------------------------------------------
  797.  
  798. function VectorReflect(const V, N: TAffineVector): TAffineVector; assembler; register;
  799.  
  800. // reflects vector V against N (assumes N is normalized)
  801. // EAX contains address of V
  802. // EDX contains address of N
  803. // ECX contains address of the result
  804.  
  805. //var Dot : Single;
  806.  
  807. asm
  808.    {Dot := VectorAffineDotProduct(V, N);
  809.    Result[X] := V[X]-2 * Dot * N[X];
  810.    Result[Y] := V[Y]-2 * Dot * N[Y];
  811.    Result[Z] := V[Z]-2 * Dot * N[Z];}
  812.    
  813.                    CALL VectorAffineDotProduct   // dot is now in ST(0)
  814.                    FCHS                          // -dot
  815.                    FADD ST, ST                   // -dot * 2
  816.                    FLD DWORD PTR [EDX]           // ST := N[X]
  817.                    FMUL ST, ST(1)                // ST := -2 * dot * N[X]
  818.                    FADD DWORD PTR[EAX]           // ST := V[X] - 2 * dot * N[X]
  819.                    FSTP DWORD PTR [ECX]          // store result
  820.                    FLD DWORD PTR [EDX + 4]       // etc.
  821.                    FMUL ST, ST(1)
  822.                    FADD DWORD PTR[EAX + 4]
  823.                    FSTP DWORD PTR [ECX + 4]
  824.                    FLD DWORD PTR [EDX + 8]
  825.                    FMUL ST, ST(1)
  826.                    FADD DWORD PTR[EAX + 8]
  827.                    FSTP DWORD PTR [ECX + 8]
  828.                    FSTP ST                       // clean FPU stack
  829. end;
  830.  
  831. // VectorRotate
  832. //
  833. procedure VectorRotate(var Vector: TVector4f; const Axis: TVector3f; Angle: Single);
  834. var
  835.    rotMatrix : TMatrix4f;
  836. begin
  837.    rotMatrix:=CreateRotationMatrix(Axis, Angle);
  838.    vector:=VectorTransform(Vector, RotMatrix);
  839. end;
  840.  
  841. // VectorScale
  842. //
  843. procedure VectorScale(V: array of Single; Factor: Single); assembler; register;
  844.  
  845. // EAX contains address of V
  846. // EDX contains highest index in V
  847. // Factor is located on the stack
  848.  
  849. asm
  850.   {for I := Low(V) to High(V) do V[I] := V[I] * Factor;}
  851.  
  852.               FLD DWORD PTR [Factor]        // load factor
  853. @@Loop:       FLD DWORD PTR [EAX + 4 * EDX] // load a component
  854.               FMUL ST, ST(1)                // multiply it with the factor
  855.               WAIT
  856.               FSTP DWORD PTR [EAX + 4 * EDX] // store the result
  857.               DEC EDX                       // do the entire array
  858.               JNS @@Loop
  859.               FSTP ST(0)                    // clean the FPU stack
  860. end;
  861.  
  862. //----------------------------------------------------------------------------------------------------------------------
  863.  
  864. procedure VectorNegate(V: array of Single); assembler; register;
  865.  
  866. // returns a negated vector
  867. // EAX contains address of V
  868. // EDX contains highest index in V
  869.  
  870. asm
  871.   {V[X] := -V[X];
  872.   V[Y] := -V[Y];
  873.   V[Z] := -V[Z];}
  874.  
  875. @@Loop:       FLD DWORD PTR [EAX + 4 * EDX]
  876.               FCHS
  877.               WAIT
  878.               FSTP DWORD PTR [EAX + 4 * EDX]
  879.               DEC EDX
  880.               JNS @@Loop
  881. end;
  882.  
  883. //----------------------------------------------------------------------------------------------------------------------
  884.  
  885. function VectorAdd(const V1, V2: TVector): TVector; register;
  886. begin
  887.   Result[X] := V1[X] + V2[X];
  888.   Result[Y] := V1[Y] + V2[Y];
  889.   Result[Z] := V1[Z] + V2[Z];
  890.   Result[W] := V1[W] + V2[W];
  891. end;
  892.  
  893. // VectorAffineAdd
  894. //
  895. function VectorAffineAdd(const V1, V2: TAffineVector): TAffineVector; register;
  896. begin
  897.    Result[X] := V1[X] + V2[X];
  898.    Result[Y] := V1[Y] + V2[Y];
  899.    Result[Z] := V1[Z] + V2[Z];
  900. end;
  901.  
  902. // VectorSubtract
  903. //
  904. function VectorSubtract(const V1, V2: TVector): TVector; register;
  905. begin
  906.   Result[X] := V1[X] - V2[X];
  907.   Result[Y] := V1[Y] - V2[Y];
  908.   Result[Z] := V1[Z] - V2[Z];
  909.   Result[W] := V1[W] - V2[W];
  910. end;
  911.  
  912. // VectorSpacing
  913. //
  914. function  VectorSpacing(const V1, V2: TVector): Single;
  915. begin
  916.    Result:=Abs(V1[X]-V2[X])+Abs(V1[Y]-V2[Y])+Abs(V1[Z]-V2[Z])+Abs(V1[W]-V2[W]);
  917. end;
  918.  
  919. // VectorEquals (hmg vector)
  920. //
  921. function  VectorEquals(const V1, V2: TVector) : Boolean;
  922. begin
  923.    Result:=((V1[0]=V2[0]) and (V1[1]=V2[1]) and (V1[2]=V2[2]) and (V1[3]=V2[3]));
  924. end;
  925.  
  926. // VectorEquals (affine vector)
  927. //
  928. function  VectorEquals(const V1, V2: TAffineVector) : Boolean;
  929. begin
  930.    Result:=((V1[0]=V2[0]) and (V1[1]=V2[1]) and (V1[2]=V2[2]));
  931. end;
  932.  
  933. //----------------------------------------------------------------------------------------------------------------------
  934.  
  935. function VectorDotProduct(const V1, V2: TVector): Single; register;
  936.  
  937. begin
  938.   Result := V1[X] * V2[X] + V1[Y] * V2[Y] + V1[Z] * V2[Z] + V1[W] * V2[W];
  939. end;
  940.  
  941. //----------------------------------------------------------------------------------------------------------------------
  942.  
  943. function VectorAffineDotProduct(const V1, V2: TAffineVector): Single; assembler; register;
  944.  
  945. // calculates the dot product between V1 and V2
  946. // EAX contains address of V1
  947. // EDX contains address of V2
  948. // result is stored in ST(0)
  949.  
  950. asm
  951.   //Result := V1[X] * V2[X] + V1[Y] * V2[Y] + V1[Z] * V2[Z];
  952.  
  953.                    FLD DWORD PTR [EAX]
  954.                    FMUL DWORD PTR [EDX]
  955.                    FLD DWORD PTR [EAX + 4]
  956.                    FMUL DWORD PTR [EDX + 4]
  957.                    FADDP
  958.                    FLD DWORD PTR [EAX + 8]
  959.                    FMUL DWORD PTR [EDX + 8]
  960.                    FADDP
  961. end;
  962.  
  963. // VectorCrossProduct
  964. //
  965. function VectorCrossProduct(const V1, V2: TAffineVector): TAffineVector;
  966.  
  967. // Temp is necessary because
  968. // either V1 or V2 could also be the result vector
  969. //
  970. // EAX contains address of V1
  971. // EDX contains address of V2
  972. // ECX contains address of result
  973.  
  974. var Temp: TAffineVector;
  975.  
  976. asm
  977.   {Temp[X] := V1[Y] * V2[Z]-V1[Z] * V2[Y];
  978.   Temp[Y] := V1[Z] * V2[X]-V1[X] * V2[Z];
  979.   Temp[Z] := V1[X] * V2[Y]-V1[Y] * V2[X];
  980.   Result := Temp;}
  981.  
  982.    PUSH EBX                      // save EBX, must be restored to original value
  983.    LEA EBX, [Temp]
  984.    FLD DWORD PTR [EDX + 8]       // first load both vectors onto FPU register stack
  985.    FLD DWORD PTR [EDX + 4]
  986.    FLD DWORD PTR [EDX + 0]
  987.    FLD DWORD PTR [EAX + 8]
  988.    FLD DWORD PTR [EAX + 4]
  989.    FLD DWORD PTR [EAX + 0]
  990.  
  991.    FLD ST(1)                     // ST(0) := V1[Y]
  992.    FMUL ST, ST(6)                // ST(0) := V1[Y] * V2[Z]
  993.    FLD ST(3)                     // ST(0) := V1[Z]
  994.    FMUL ST, ST(6)                // ST(0) := V1[Z] * V2[Y]
  995.    FSUBP ST(1), ST               // ST(0) := ST(1)-ST(0)
  996.    FSTP DWORD [EBX]              // Temp[X] := ST(0)
  997.    FLD ST(2)                     // ST(0) := V1[Z]
  998.    FMUL ST, ST(4)                // ST(0) := V1[Z] * V2[X]
  999.    FLD ST(1)                     // ST(0) := V1[X]
  1000.    FMUL ST, ST(7)                // ST(0) := V1[X] * V2[Z]
  1001.    FSUBP ST(1), ST               // ST(0) := ST(1)-ST(0)
  1002.    FSTP DWORD [EBX + 4]          // Temp[Y] := ST(0)
  1003.    FLD ST                        // ST(0) := V1[X]
  1004.    FMUL ST, ST(5)                // ST(0) := V1[X] * V2[Y]
  1005.    FLD ST(2)                     // ST(0) := V1[Y]
  1006.    FMUL ST, ST(5)                // ST(0) := V1[Y] * V2[X]
  1007.    FSUBP ST(1), ST               // ST(0) := ST(1)-ST(0)
  1008.    FSTP DWORD [EBX + 8]          // Temp[Z] := ST(0)
  1009.    FSTP ST(0)                    // clear FPU register stack
  1010.    FSTP ST(0)
  1011.    FSTP ST(0)
  1012.    FSTP ST(0)
  1013.    FSTP ST(0)
  1014.    FSTP ST(0)
  1015.    MOV EAX, [EBX]                // copy Temp to Result
  1016.    MOV [ECX], EAX
  1017.    MOV EAX, [EBX + 4]
  1018.    MOV [ECX + 4], EAX
  1019.    MOV EAX, [EBX + 8]
  1020.    MOV [ECX + 8], EAX
  1021.    POP EBX
  1022. end;
  1023.  
  1024. // VectorCrossProduct
  1025. //
  1026. function  VectorCrossProduct(const V1, V2: TVector): TVector;
  1027. begin
  1028.    Result[X] := V1[Y] * V2[Z]-V1[Z] * V2[Y];
  1029.    Result[Y] := V1[Z] * V2[X]-V1[X] * V2[Z];
  1030.    Result[Z] := V1[X] * V2[Y]-V1[Y] * V2[X];
  1031.    Result[W] := 0;
  1032. end;
  1033.  
  1034. // TAffineVector
  1035. //
  1036. function VectorPerpendicular(const V, N: TAffineVector): TAffineVector;
  1037. var
  1038.    dot : Single;
  1039. begin
  1040.    dot := VectorAffineDotProduct(V, N);
  1041.    Result[X] := V[X]-Dot * N[X];
  1042.    Result[Y] := V[Y]-Dot * N[Y];
  1043.    Result[Z] := V[Z]-Dot * N[Z];
  1044. end;
  1045.  
  1046. // VectorTransform
  1047. //
  1048. function VectorTransform(const V: TVector4f; const M: TMatrix): TVector4f; register;
  1049. var
  1050.    TV: TVector4f;
  1051. begin
  1052.    TV[X] := V[X] * M[X, X] + V[Y] * M[Y, X] + V[Z] * M[Z, X] + V[W] * M[W, X];
  1053.    TV[Y] := V[X] * M[X, Y] + V[Y] * M[Y, Y] + V[Z] * M[Z, Y] + V[W] * M[W, Y];
  1054.    TV[Z] := V[X] * M[X, Z] + V[Y] * M[Y, Z] + V[Z] * M[Z, Z] + V[W] * M[W, Z];
  1055.    TV[W] := V[X] * M[X, W] + V[Y] * M[Y, W] + V[Z] * M[Z, W] + V[W] * M[W, W];
  1056.    Result := TV
  1057. end;
  1058.  
  1059. // VectorTransform
  1060. //
  1061. function  VectorTransform(const V: TVector3f; const M: TMatrix): TVector3f;
  1062. var
  1063.    TV: TVector3f;
  1064. begin
  1065.    TV[X] := V[X] * M[X, X] + V[Y] * M[Y, X] + V[Z] * M[Z, X] + M[W, X];
  1066.    TV[Y] := V[X] * M[X, Y] + V[Y] * M[Y, Y] + V[Z] * M[Z, Y] + M[W, Y];
  1067.    TV[Z] := V[X] * M[X, Z] + V[Y] * M[Y, Z] + V[Z] * M[Z, Z] + M[W, Z];
  1068.    Result := TV;
  1069. end;
  1070.  
  1071. // VectorAffineTransform
  1072. //
  1073. function VectorAffineTransform(const V: TAffineVector; const M: TAffineMatrix): TAffineVector; register;
  1074. var
  1075.    TV: TAffineVector;
  1076. begin
  1077.    TV[X] := V[X] * M[X, X] + V[Y] * M[Y, X] + V[Z] * M[Z, X];
  1078.    TV[Y] := V[X] * M[X, Y] + V[Y] * M[Y, Y] + V[Z] * M[Z, Y];
  1079.    TV[Z] := V[X] * M[X, Z] + V[Y] * M[Y, Z] + V[Z] * M[Z, Z];
  1080.    Result := TV;
  1081. end;
  1082.  
  1083. //----------------------------------------------------------------------------------------------------------------------
  1084.  
  1085. function PointInPolygon(xp, yp : array of Single; x, y: Single) : Boolean;
  1086. // The code below is from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
  1087. // with some minor modifications for speed.  It returns 1 for strictly
  1088. // interior points, 0 for strictly exterior, and 0 or 1 for points on
  1089. // the boundary.
  1090. // This code is not yet tested!
  1091. var
  1092.    i, j: Integer;
  1093. begin
  1094.    Result:=False;
  1095.    if High(XP)=High(YP) then begin
  1096.       j:=High(XP);
  1097.       for i:=0 to High(XP) do begin
  1098.          if ( ( ((yp[I]<=y) and (y<yp[J])) or ((yp[J]<=y) and (y<yp[I])) )
  1099.              and (x<(xp[J]-xp[I])*(y-yp[I])/(yp[J]-yp[I])+xp[I])) then
  1100.             Result := not Result;
  1101.          j:=i;
  1102.          //jh mod
  1103.       end;
  1104.    end;
  1105. end;
  1106.  
  1107. //----------------------------------------------------------------------------------------------------------------------
  1108.  
  1109. function QuaternionConjugate(const Q: TQuaternion): TQuaternion; assembler;
  1110.  
  1111. // returns the conjugate of a quaternion
  1112. // EAX contains address of Q
  1113. // EDX contains address of result
  1114.  
  1115. asm
  1116.               FLD DWORD PTR [EAX]
  1117.               FCHS
  1118.               WAIT
  1119.               FSTP DWORD PTR [EDX]
  1120.               FLD DWORD PTR [EAX + 4]
  1121.               FCHS
  1122.               WAIT
  1123.               FSTP DWORD PTR [EDX + 4]
  1124.               FLD DWORD PTR [EAX + 8]
  1125.               FCHS
  1126.               WAIT
  1127.               FSTP DWORD PTR [EDX + 8]
  1128.               MOV EAX, [EAX + 12]
  1129.               MOV [EDX + 12], EAX
  1130. end;
  1131.  
  1132. //----------------------------------------------------------------------------------------------------------------------
  1133.  
  1134. function QuaternionFromPoints(const V1, V2: TAffineVector): TQuaternion; assembler;
  1135.  
  1136. // constructs a unit quaternion from two points on unit sphere
  1137. // EAX contains address of V1
  1138. // ECX contains address to result
  1139. // EDX contains address of V2
  1140.  
  1141. asm
  1142.   {Result.ImagPart := VectorCrossProduct(V1, V2);
  1143.    Result.RealPart :=  Sqrt((VectorAffineDotProduct(V1, V2) + 1)/2);}
  1144.  
  1145.               PUSH EAX
  1146.               CALL VectorCrossProduct       // determine axis to rotate about
  1147.               POP EAX
  1148.               FLD1                          // prepare next calculation
  1149.               Call VectorAffineDotProduct   // calculate cos(angle between V1 and V2)
  1150.               FADD ST, ST(1)                // transform angle to angle/2 by: cos(a/2)=sqrt((1 + cos(a))/2)
  1151.               FXCH ST(1)
  1152.               FADD ST, ST
  1153.               FDIVP ST(1), ST
  1154.               FSQRT
  1155.               FSTP DWORD PTR [ECX + 12]     // Result.RealPart := ST(0)
  1156. end;
  1157.  
  1158. //----------------------------------------------------------------------------------------------------------------------
  1159.  
  1160. function QuaternionMultiply(const qL, qR: TQuaternion): TQuaternion;
  1161.  
  1162. // Returns quaternion product qL * qR.  Note: order is important!
  1163. // To combine rotations, use the product QuaternionMuliply(qSecond, qFirst),
  1164. // which gives the effect of rotating by qFirst then qSecond.
  1165.  
  1166. var Temp : TQuaternion;
  1167.  
  1168. begin
  1169.   Temp.RealPart := qL.RealPart * qR.RealPart - qL.ImagPart[X] * qR.ImagPart[X] -
  1170.                    qL.ImagPart[Y] * qR.ImagPart[Y] - qL.ImagPart[Z] * qR.ImagPart[Z];
  1171.   Temp.ImagPart[X] := qL.RealPart * qR.ImagPart[X] + qL.ImagPart[X] * qR.RealPart +
  1172.                       qL.ImagPart[Y] * qR.ImagPart[Z] - qL.ImagPart[Z] * qR.ImagPart[Y];
  1173.   Temp.ImagPart[Y] := qL.RealPart * qR.ImagPart[Y] + qL.ImagPart[Y] * qR.RealPart +
  1174.                       qL.ImagPart[Z] * qR.ImagPart[X] - qL.ImagPart[X] * qR.ImagPart[Z];
  1175.   Temp.ImagPart[Z] := qL.RealPart * qR.ImagPart[Z] + qL.ImagPart[Z] * qR.RealPart +
  1176.                       qL.ImagPart[X] * qR.ImagPart[Y] - qL.ImagPart[Y] * qR.ImagPart[X];
  1177.   Result := Temp;
  1178. end;
  1179.  
  1180. //----------------------------------------------------------------------------------------------------------------------
  1181.  
  1182. function QuaternionToMatrix(const Q: TQuaternion): TMatrix;
  1183.  
  1184. // Constructs rotation matrix from (possibly non-unit) quaternion.
  1185. // Assumes matrix is used to multiply column vector on the left:
  1186. // vnew = mat vold.  Works correctly for right-handed coordinate system
  1187. // and right-handed rotations.
  1188.  
  1189. // Essentially, this function is the same as CreateRotationMatrix and you can consider it as
  1190. // being for reference here.
  1191.  
  1192. {var Norm, S,
  1193.     XS, YS, ZS,
  1194.     WX, WY, WZ,
  1195.     XX, XY, XZ,
  1196.     YY, YZ, ZZ   : Single;
  1197.  
  1198. begin
  1199.   Norm := Q.Vector[X] * Q.Vector[X] + Q.Vector[Y] * Q.Vector[Y] + Q.Vector[Z] * Q.Vector[Z] + Q.RealPart * Q.RealPart;
  1200.   if Norm > 0 then S := 2 / Norm
  1201.               else S := 0;
  1202.  
  1203.   XS := Q.Vector[X] * S;   YS := Q.Vector[Y] * S;   ZS := Q.Vector[Z] * S;
  1204.   WX := Q.RealPart * XS;   WY := Q.RealPart * YS;   WZ := Q.RealPart * ZS;
  1205.   XX := Q.Vector[X] * XS;  XY := Q.Vector[X] * YS;  XZ := Q.Vector[X] * ZS;
  1206.   YY := Q.Vector[Y] * YS;  YZ := Q.Vector[Y] * ZS;  ZZ := Q.Vector[Z] * ZS;
  1207.  
  1208.   Result[X, X] := 1 - (YY + ZZ); Result[Y, X] := XY + WZ;       Result[Z, X] := XZ - WY;       Result[W, X] := 0;
  1209.   Result[X, Y] := XY - WZ;       Result[Y, Y] := 1 - (XX + ZZ); Result[Z, Y] := YZ + WX;       Result[W, Y] := 0;
  1210.   Result[X, Z] := XZ + WY;       Result[Y, Z] := YZ - WX;       Result[Z, Z] := 1 - (XX + YY); Result[W, Z] := 0;
  1211.   Result[X, W] := 0;             Result[Y, W] := 0;             Result[Z, W] := 0;             Result[W, W] := 1;}
  1212.  
  1213. var
  1214.   V: TAffineVector;
  1215.   SinA, CosA,
  1216.   A, B, C: Extended;
  1217.  
  1218. begin
  1219.   V := Q.ImagPart;
  1220.   VectorNormalize(V);
  1221.   SinCos(Q.RealPart / 2, SinA, CosA);
  1222.   A := V[X] * SinA;
  1223.   B := V[Y] * SinA;
  1224.   C := V[Z] * SinA;
  1225.  
  1226.   Result := IdentityMatrix;
  1227.   Result[X, X] := 1 - 2 * B * B - 2 * C * C;
  1228.   Result[X, Y] := 2 * A * B - 2 * CosA * C;
  1229.   Result[X, Z] := 2 * A * C + 2 * CosA * B;
  1230.  
  1231.   Result[Y, X] := 2 * A * B + 2 * CosA * C;
  1232.   Result[Y, Y] := 1 - 2 * A * A - 2 * C * C;
  1233.   Result[Y, Z] := 2 * B * C - 2 * CosA * A;
  1234.  
  1235.   Result[Z, X] := 2 * A * C - 2 * CosA * B;
  1236.   Result[Z, Y] := 2 * B * C + 2 * CosA * A;
  1237.   Result[Z, Z] := 1 - 2 * A * A - 2 * B * B;
  1238. end;
  1239.  
  1240. //----------------------------------------------------------------------------------------------------------------------
  1241.  
  1242. procedure QuaternionToPoints(const Q: TQuaternion; var ArcFrom, ArcTo: TAffineVector); register;
  1243.  
  1244. // converts a unit quaternion into two points on a unit sphere
  1245.  
  1246. var S: Single;
  1247.  
  1248. begin
  1249.   S := Sqrt(Q.ImagPart[X] * Q.ImagPart[X] + Q.ImagPart[Y] * Q.ImagPart[Y]);
  1250.   if S = 0 then ArcFrom := MakeAffineVector([0, 1, 0])
  1251.            else ArcFrom := MakeAffineVector([-Q.ImagPart[Y] / S, Q.ImagPart[X] / S, 0]);
  1252.   ArcTo[X] := Q.RealPart * ArcFrom[X] - Q.ImagPart[Z] * ArcFrom[Y];
  1253.   ArcTo[Y] := Q.RealPart * ArcFrom[Y] + Q.ImagPart[Z] * ArcFrom[X];
  1254.   ArcTo[Z] := Q.ImagPart[X] * ArcFrom[Y] - Q.ImagPart[Y] * ArcFrom[X];
  1255.   if Q.RealPart < 0 then ArcFrom := MakeAffineVector([-ArcFrom[X], -ArcFrom[Y], 0]);
  1256. end;
  1257.  
  1258. //----------------------------------------------------------------------------------------------------------------------
  1259.  
  1260. function MatrixAffineDeterminant(const M: TAffineMatrix): Single; register;
  1261.  
  1262. // determinant of a 3x3 matrix
  1263.  
  1264. begin
  1265.   Result := M[X, X] * (M[Y, Y] * M[Z, Z] - M[Z, Y] * M[Y, Z]) -
  1266.             M[X, Y] * (M[Y, X] * M[Z, Z] - M[Z, X] * M[Y, Z]) +
  1267.             M[X, Z] * (M[Y, X] * M[Z, Y] - M[Z, X] * M[Y, Y]);
  1268. end;
  1269.  
  1270. //----------------------------------------------------------------------------------------------------------------------
  1271.  
  1272. function MatrixDetInternal(a1, a2, a3, b1, b2, b3, c1, c2, c3: Single): Single;
  1273.  
  1274. // internal version for the determinant of a 3x3 matrix
  1275.  
  1276. begin
  1277.   Result := a1 * (b2 * c3 - b3 * c2) -
  1278.             b1 * (a2 * c3 - a3 * c2) +
  1279.             c1 * (a2 * b3 - a3 * b2);
  1280. end;
  1281.  
  1282. //----------------------------------------------------------------------------------------------------------------------
  1283.  
  1284. procedure MatrixAdjoint(var M: TMatrix); register;
  1285.  
  1286. // Adjoint of a 4x4 matrix - used in the computation of the inverse
  1287. // of a 4x4 matrix
  1288.  
  1289. var a1, a2, a3, a4, 
  1290.     b1, b2, b3, b4, 
  1291.     c1, c2, c3, c4, 
  1292.     d1, d2, d3, d4: Single;
  1293.  
  1294.  
  1295. begin
  1296.     a1 :=  M[X, X]; b1 :=  M[X, Y];
  1297.     c1 :=  M[X, Z]; d1 :=  M[X, W];
  1298.     a2 :=  M[Y, X]; b2 :=  M[Y, Y];
  1299.     c2 :=  M[Y, Z]; d2 :=  M[Y, W];
  1300.     a3 :=  M[Z, X]; b3 :=  M[Z, Y];
  1301.     c3 :=  M[Z, Z]; d3 :=  M[Z, W];
  1302.     a4 :=  M[W, X]; b4 :=  M[W, Y];
  1303.     c4 :=  M[W, Z]; d4 :=  M[W, W];
  1304.  
  1305.     // row column labeling reversed since we transpose rows & columns
  1306.     M[X, X] :=  MatrixDetInternal(b2, b3, b4, c2, c3, c4, d2, d3, d4);
  1307.     M[Y, X] := -MatrixDetInternal(a2, a3, a4, c2, c3, c4, d2, d3, d4);
  1308.     M[Z, X] :=  MatrixDetInternal(a2, a3, a4, b2, b3, b4, d2, d3, d4);
  1309.     M[W, X] := -MatrixDetInternal(a2, a3, a4, b2, b3, b4, c2, c3, c4);
  1310.  
  1311.     M[X, Y] := -MatrixDetInternal(b1, b3, b4, c1, c3, c4, d1, d3, d4);
  1312.     M[Y, Y] :=  MatrixDetInternal(a1, a3, a4, c1, c3, c4, d1, d3, d4);
  1313.     M[Z, Y] := -MatrixDetInternal(a1, a3, a4, b1, b3, b4, d1, d3, d4);
  1314.     M[W, Y] :=  MatrixDetInternal(a1, a3, a4, b1, b3, b4, c1, c3, c4);
  1315.  
  1316.     M[X, Z] :=  MatrixDetInternal(b1, b2, b4, c1, c2, c4, d1, d2, d4);
  1317.     M[Y, Z] := -MatrixDetInternal(a1, a2, a4, c1, c2, c4, d1, d2, d4);
  1318.     M[Z, Z] :=  MatrixDetInternal(a1, a2, a4, b1, b2, b4, d1, d2, d4);
  1319.     M[W, Z] := -MatrixDetInternal(a1, a2, a4, b1, b2, b4, c1, c2, c4);
  1320.  
  1321.     M[X, W] := -MatrixDetInternal(b1, b2, b3, c1, c2, c3, d1, d2, d3);
  1322.     M[Y, W] :=  MatrixDetInternal(a1, a2, a3, c1, c2, c3, d1, d2, d3);
  1323.     M[Z, W] := -MatrixDetInternal(a1, a2, a3, b1, b2, b3, d1, d2, d3);
  1324.     M[W, W] :=  MatrixDetInternal(a1, a2, a3, b1, b2, b3, c1, c2, c3);
  1325. end;
  1326.  
  1327. //----------------------------------------------------------------------------------------------------------------------
  1328.  
  1329. function MatrixDeterminant(const M: TMatrix): Single; register;
  1330.  
  1331. // Determinant of a 4x4 matrix
  1332.  
  1333. var a1, a2, a3, a4, 
  1334.     b1, b2, b3, b4,
  1335.     c1, c2, c3, c4, 
  1336.     d1, d2, d3, d4  : Single;
  1337.  
  1338. begin
  1339.   a1 := M[X, X];  b1 := M[X, Y];  c1 := M[X, Z];  d1 := M[X, W];
  1340.   a2 := M[Y, X];  b2 := M[Y, Y];  c2 := M[Y, Z];  d2 := M[Y, W];
  1341.   a3 := M[Z, X];  b3 := M[Z, Y];  c3 := M[Z, Z];  d3 := M[Z, W];
  1342.   a4 := M[W, X];  b4 := M[W, Y];  c4 := M[W, Z];  d4 := M[W, W];
  1343.  
  1344.   Result := a1 * MatrixDetInternal(b2, b3, b4, c2, c3, c4, d2, d3, d4) -
  1345.             b1 * MatrixDetInternal(a2, a3, a4, c2, c3, c4, d2, d3, d4) +
  1346.             c1 * MatrixDetInternal(a2, a3, a4, b2, b3, b4, d2, d3, d4) -
  1347.             d1 * MatrixDetInternal(a2, a3, a4, b2, b3, b4, c2, c3, c4);
  1348. end;
  1349.  
  1350. //----------------------------------------------------------------------------------------------------------------------
  1351.  
  1352. procedure MatrixScale(var M: TMatrix; Factor: Single); register;
  1353.  
  1354. // multiplies all elements of a 4x4 matrix with a factor
  1355.  
  1356. var I, J: Integer;
  1357.  
  1358. begin
  1359.   for I := 0 to 3 do
  1360.     for J := 0 to 3 do M[I, J] := M[I, J] * Factor;
  1361. end;
  1362.  
  1363. //----------------------------------------------------------------------------------------------------------------------
  1364.  
  1365. procedure MatrixInvert(var M: TMatrix); register;
  1366.  
  1367. // finds the inverse of a 4x4 matrix
  1368.  
  1369. var Det: Single;
  1370.  
  1371. begin
  1372.   Det := MatrixDeterminant(M);
  1373.   if Abs(Det) < EPSILON then M := IdentityMatrix
  1374.                         else
  1375.   begin
  1376.     MatrixAdjoint(M);
  1377.     MatrixScale(M, 1 / Det);
  1378.   end;
  1379. end;
  1380.  
  1381. //----------------------------------------------------------------------------------------------------------------------
  1382.  
  1383. procedure MatrixTranspose(var M: TMatrix); register;
  1384.  
  1385. // computes transpose of 4x4 matrix
  1386.  
  1387. var I, J: Integer;
  1388.     TM: TMatrix;
  1389.  
  1390. begin
  1391.   for I := 0 to 3 do
  1392.     for J := 0 to 3 do TM[J, I] := M[I, J];
  1393.   M := TM;
  1394. end;
  1395.  
  1396. //----------------------------------------------------------------------------------------------------------------------
  1397.  
  1398. procedure MatrixAffineTranspose(var M: TAffineMatrix); register;
  1399.  
  1400. // computes transpose of 3x3 matrix
  1401.  
  1402. var I, J: Integer;
  1403.     TM: TAffineMatrix;
  1404.  
  1405. begin
  1406.   for I := 0 to 2 do
  1407.     for J := 0 to 2 do TM[J, I] := M[I, J];
  1408.   M := TM;
  1409. end;
  1410.  
  1411. //----------------------------------------------------------------------------------------------------------------------
  1412.  
  1413. function MatrixMultiply(const M1, M2: TMatrix): TMatrix; register;
  1414.  
  1415. // multiplies two 4x4 matrices
  1416.  
  1417. var I, J: Integer;
  1418.     TM: TMatrix;
  1419.  
  1420. begin
  1421.   for I := 0 to 3 do
  1422.     for J := 0 to 3 do
  1423.       TM[I, J] := M1[I, X] * M2[X, J] +
  1424.                   M1[I, Y] * M2[Y, J] +
  1425.                   M1[I, Z] * M2[Z, J] +
  1426.                   M1[I, W] * M2[W, J];
  1427.   Result := TM;
  1428. end;
  1429.  
  1430. //----------------------------------------------------------------------------------------------------------------------
  1431.  
  1432. function CreateRotationMatrix(const Axis: TVector3f; Angle: Single): TMatrix; register;
  1433.  
  1434. // Creates a rotation matrix along the given Axis by the given Angle in radians.
  1435.  
  1436. var cosine,
  1437.     sine,
  1438.     Len,
  1439.     one_minus_cosine: Extended;
  1440.  
  1441. begin
  1442.   SinCos(Angle, Sine, Cosine);
  1443.   one_minus_cosine := 1 - cosine;
  1444.   Len := VectorNormalize(Axis);
  1445.  
  1446.   if Len = 0 then Result := IdentityMatrix
  1447.              else
  1448.   begin
  1449.     Result[X, X] := (one_minus_cosine * Sqr(Axis[0])) + Cosine;
  1450.     Result[X, Y] := (one_minus_cosine * Axis[0] * Axis[1]) - (Axis[2] * Sine);
  1451.     Result[X, Z] := (one_minus_cosine * Axis[2] * Axis[0]) + (Axis[1] * Sine);
  1452.     Result[X, W] := 0;
  1453.  
  1454.     Result[Y, X] := (one_minus_cosine * Axis[0] * Axis[1]) + (Axis[2] * Sine);
  1455.     Result[Y, Y] := (one_minus_cosine * Sqr(Axis[1])) + Cosine;
  1456.     Result[Y, Z] := (one_minus_cosine * Axis[1] * Axis[2]) - (Axis[0] * Sine);
  1457.     Result[Y, W] := 0;
  1458.  
  1459.     Result[Z, X] := (one_minus_cosine * Axis[2] * Axis[0]) - (Axis[1] * Sine);
  1460.     Result[Z, Y] := (one_minus_cosine * Axis[1] * Axis[2]) + (Axis[0] * Sine);
  1461.     Result[Z, Z] := (one_minus_cosine * Sqr(Axis[2])) + Cosine;
  1462.     Result[Z, W] := 0;
  1463.  
  1464.     Result[W, X] := 0;
  1465.     Result[W, Y] := 0;
  1466.     Result[W, Z] := 0;
  1467.     Result[W, W] := 1;
  1468.   end;
  1469. end;
  1470.  
  1471. //----------------------------------------------------------------------------------------------------------------------
  1472.  
  1473. function ConvertRotation(const Angles: TAffineVector): TVector; register;
  1474.  
  1475. { Turn a triplet of rotations about x, y, and z (in that order) into an
  1476.    equivalent rotation around a single axis (all in radians).
  1477.  
  1478.    Rotation of the Angle t about the axis (X, Y, Z) is given by:
  1479.  
  1480.      | X^2 + (1-X^2) Cos(t),    XY(1-Cos(t))  +  Z Sin(t), XZ(1-Cos(t))-Y Sin(t) |
  1481.  M = | XY(1-Cos(t))-Z Sin(t), Y^2 + (1-Y^2) Cos(t),      YZ(1-Cos(t)) + X Sin(t) |
  1482.      | XZ(1-Cos(t)) + Y Sin(t), YZ(1-Cos(t))-X Sin(t),   Z^2 + (1-Z^2) Cos(t)    |
  1483.  
  1484.    Rotation about the three axes (Angles a1, a2, a3) can be represented as
  1485.    the product of the individual rotation matrices:
  1486.  
  1487.       | 1  0       0       | | Cos(a2) 0 -Sin(a2) | |  Cos(a3) Sin(a3) 0 |
  1488.       | 0  Cos(a1) Sin(a1) | * | 0       1  0       | * | -Sin(a3) Cos(a3) 0 |
  1489.       | 0 -Sin(a1) Cos(a1) | | Sin(a2) 0  Cos(a2) | |  0       0       1 |
  1490.          Mx                       My                     Mz
  1491.  
  1492.    We now want to solve for X, Y, Z, and t given 9 equations in 4 unknowns.
  1493.    Using the diagonal elements of the two matrices, we get:
  1494.  
  1495.       X^2 + (1-X^2) Cos(t) = M[0][0]
  1496.       Y^2 + (1-Y^2) Cos(t) = M[1][1]
  1497.       Z^2 + (1-Z^2) Cos(t) = M[2][2]
  1498.  
  1499.    Adding the three equations, we get:
  1500.  
  1501.       X^2  +  Y^2  +  Z^2 - (M[0][0]  +  M[1][1]  +  M[2][2]) =
  1502.      - (3 - X^2 - Y^2 - Z^2) Cos(t)
  1503.  
  1504.    Since (X^2  +  Y^2  +  Z^2) = 1, we can rewrite as:
  1505.  
  1506.       Cos(t) = (1 - (M[0][0]  +  M[1][1]  +  M[2][2])) / 2
  1507.  
  1508.    Solving for t, we get:
  1509.  
  1510.       t = Acos(((M[0][0]  +  M[1][1]  +  M[2][2]) - 1) / 2)
  1511.  
  1512.     We can substitute t into the equations for X^2, Y^2, and Z^2 above
  1513.     to get the values for X, Y, and Z.  To find the proper signs we note
  1514.     that:
  1515.  
  1516.     2 X Sin(t) = M[1][2] - M[2][1]
  1517.     2 Y Sin(t) = M[2][0] - M[0][2]
  1518.     2 Z Sin(t) = M[0][1] - M[1][0]
  1519. }
  1520.  
  1521. var Axis1, Axis2: TVector3f;
  1522.     M, M1, M2: TMatrix;
  1523.     cost, cost1, 
  1524.     sint, 
  1525.     s1, s2, s3: Single;
  1526.     I: Integer;
  1527.  
  1528.  
  1529. begin
  1530.    // see if we are only rotating about a single Axis
  1531.    if Abs(Angles[X]) < EPSILON then begin
  1532.       if Abs(Angles[Y]) < EPSILON then begin
  1533.          Result := MakeVector([0, 0, 1, Angles[Z]]);
  1534.          Exit;
  1535.       end else if Abs(Angles[Z]) < EPSILON then begin
  1536.          Result := MakeVector([0, 1, 0, Angles[Y]]);
  1537.          Exit;
  1538.       end
  1539.    end else if (Abs(Angles[Y]) < EPSILON) and (Abs(Angles[Z]) < EPSILON) then begin
  1540.       Result := MakeVector([1, 0, 0, Angles[X]]);
  1541.       Exit;
  1542.    end;
  1543.  
  1544.    // make the rotation matrix
  1545.    Axis1 := XVector;
  1546.    M := CreateRotationMatrix(Axis1, Angles[X]);
  1547.  
  1548.    Axis2 := YVector;
  1549.    M2 := CreateRotationMatrix(Axis2, Angles[Y]);
  1550.    M1 := MatrixMultiply(M, M2);
  1551.  
  1552.    Axis2 := ZVector;
  1553.    M2 := CreateRotationMatrix(Axis2, Angles[Z]);
  1554.    M := MatrixMultiply(M1, M2);
  1555.  
  1556.    cost := ((M[X, X] + M[Y, Y] + M[Z, Z])-1) / 2;
  1557.    if cost < -1 then
  1558.       cost := -1
  1559.    else if cost > 1 - EPSILON then begin
  1560.       // Bad Angle - this would cause a crash
  1561.       Result := MakeVector([1, 0, 0, 0]);
  1562.       Exit;
  1563.    end;
  1564.  
  1565.   cost1 := 1 - cost;
  1566.   Result := Makevector([Sqrt((M[X, X]-cost) / cost1),
  1567.                       Sqrt((M[Y, Y]-cost) / cost1),
  1568.                       sqrt((M[Z, Z]-cost) / cost1),
  1569.                       arccos(cost)]);
  1570.  
  1571.   sint := 2 * Sqrt(1 - cost * cost); // This is actually 2 Sin(t)
  1572.  
  1573.   // Determine the proper signs
  1574.   for I := 0 to 7 do
  1575.   begin
  1576.     if (I and 1) > 1 then s1 := -1 else s1 := 1;
  1577.     if (I and 2) > 1 then s2 := -1 else s2 := 1;
  1578.     if (I and 4) > 1 then s3 := -1 else s3 := 1;
  1579.     if (Abs(s1 * Result[X] * sint-M[Y, Z] + M[Z, Y]) < EPSILON2) and
  1580.        (Abs(s2 * Result[Y] * sint-M[Z, X] + M[X, Z]) < EPSILON2) and
  1581.        (Abs(s3 * Result[Z] * sint-M[X, Y] + M[Y, X]) < EPSILON2) then
  1582.         begin
  1583.           // We found the right combination of signs
  1584.           Result[X] := Result[X] * s1;
  1585.           Result[Y] := Result[Y] * s2;
  1586.           Result[Z] := Result[Z] * s3;
  1587.           Exit;
  1588.         end;
  1589.   end;
  1590. end;
  1591.  
  1592. function CreateRotationMatrixX(Sine, Cosine: Single): TMatrix; register;
  1593. begin
  1594.    Result := EmptyMatrix;
  1595.    Result[X, X] := 1;
  1596.    Result[Y, Y] := Cosine;
  1597.    Result[Y, Z] := Sine;
  1598.    Result[Z, Y] := -Sine;
  1599.    Result[Z, Z] := Cosine;
  1600.    Result[W, W] := 1;
  1601. end;
  1602.  
  1603. function CreateRotationMatrixY(Sine, Cosine: Single): TMatrix; register;
  1604. begin
  1605.   Result := EmptyMatrix;
  1606.   Result[X, X] := Cosine;
  1607.   Result[X, Z] := -Sine;
  1608.   Result[Y, Y] := 1;
  1609.   Result[Z, X] := Sine;
  1610.   Result[Z, Z] := Cosine;
  1611.   Result[W, W] := 1;
  1612. end;
  1613.  
  1614. //----------------------------------------------------------------------------------------------------------------------
  1615.  
  1616. function CreateRotationMatrixZ(Sine, Cosine: Single): TMatrix; register;
  1617.  
  1618. // creates matrix for rotation about z-axis
  1619.  
  1620. begin
  1621.   Result := EmptyMatrix;
  1622.   Result[X, X] := Cosine;
  1623.   Result[X, Y] := Sine;
  1624.   Result[Y, X] := -Sine;
  1625.   Result[Y, Y] := Cosine;
  1626.   Result[Z, Z] := 1;
  1627.   Result[W, W] := 1;
  1628. end;
  1629.  
  1630. //----------------------------------------------------------------------------------------------------------------------
  1631.  
  1632. function CreateScaleMatrix(const V: TAffineVector): TMatrix; register;
  1633.  
  1634. // creates scaling matrix
  1635.  
  1636. begin
  1637.   Result := IdentityMatrix;
  1638.   Result[X, X] := V[X];
  1639.   Result[Y, Y] := V[Y];
  1640.   Result[Z, Z] := V[Z];
  1641. end;
  1642.  
  1643. //----------------------------------------------------------------------------------------------------------------------
  1644.  
  1645. function CreateTranslationMatrix(const V: TVector): TMatrix; register;
  1646.  
  1647. // creates translation matrix
  1648.  
  1649. begin
  1650.   Result := IdentityMatrix;
  1651.   Result[W, X] := V[X];
  1652.   Result[W, Y] := V[Y];
  1653.   Result[W, Z] := V[Z];
  1654.   Result[W, W] := V[W];
  1655. end;
  1656.  
  1657. //----------------------------------------------------------------------------------------------------------------------
  1658.  
  1659. function Lerp(Start, Stop, t: Single): Single;
  1660.  
  1661. // calculates linear interpolation between start and stop at point t
  1662.  
  1663. begin
  1664.   Result := Start + (Stop - Start) * t;
  1665. end;
  1666.  
  1667. //----------------------------------------------------------------------------------------------------------------------
  1668.  
  1669. function VectorAffineLerp(const V1, V2: TAffineVector; t: Single): TAffineVector;
  1670.  
  1671. // calculates linear interpolation between vector1 and vector2 at point t
  1672.  
  1673. begin
  1674.   Result[X] := Lerp(V1[X], V2[X], t);
  1675.   Result[Y] := Lerp(V1[Y], V2[Y], t);
  1676.   Result[Z] := Lerp(V1[Z], V2[Z], t);
  1677. end;
  1678.  
  1679. //----------------------------------------------------------------------------------------------------------------------
  1680.  
  1681. function VectorLerp(const V1, V2: TVector; t: Single): TVector;
  1682.  
  1683. // calculates linear interpolation between vector1 and vector2 at point t
  1684.  
  1685. begin
  1686.   Result[X] := Lerp(V1[X], V2[X], t);
  1687.   Result[Y] := Lerp(V1[Y], V2[Y], t);
  1688.   Result[Z] := Lerp(V1[Z], V2[Z], t);
  1689.   Result[W] := Lerp(V1[W], V2[W], t);
  1690. end;
  1691.  
  1692. //----------------------------------------------------------------------------------------------------------------------
  1693.  
  1694. function QuaternionSlerp(const QStart, QEnd: TQuaternion; Spin: Integer; t: Single): TQuaternion;
  1695.  
  1696. // spherical linear interpolation of unit quaternions with spins
  1697. // QStart, QEnd - start and end unit quaternions
  1698. // t            - interpolation parameter (0 to 1)
  1699. // Spin         - number of extra spin rotations to involve
  1700.  
  1701. var beta,                   // complementary interp parameter
  1702.     theta,                  // Angle between A and B
  1703.     sint, cost,             // sine, cosine of theta
  1704.     phi: Single;            // theta plus spins
  1705.     bflip: Boolean;         // use negativ t?
  1706.  
  1707.  
  1708. begin
  1709.   // cosine theta
  1710.   cost := VectorAngle(QStart.ImagPart, QEnd.ImagPart);
  1711.  
  1712.   // if QEnd is on opposite hemisphere from QStart, use -QEnd instead
  1713.   if cost < 0 then
  1714.   begin
  1715.     cost := -cost;
  1716.     bflip := True;
  1717.   end
  1718.   else bflip := False;
  1719.  
  1720.   // if QEnd is (within precision limits) the same as QStart, 
  1721.   // just linear interpolate between QStart and QEnd.
  1722.   // Can't do spins, since we don't know what direction to spin.
  1723.  
  1724.   if (1 - cost) < EPSILON then beta := 1 - t
  1725.                           else
  1726.   begin
  1727.     // normal case
  1728.     theta := arccos(cost);
  1729.     phi := theta + Spin * Pi;
  1730.     sint := sin(theta);
  1731.     beta := sin(theta - t * phi) / sint;
  1732.     t := sin(t * phi) / sint;
  1733.   end;
  1734.  
  1735.   if bflip then t := -t;
  1736.  
  1737.   // interpolate
  1738.   Result.ImagPart[X] := beta * QStart.ImagPart[X] + t * QEnd.ImagPart[X];
  1739.   Result.ImagPart[Y] := beta * QStart.ImagPart[Y] + t * QEnd.ImagPart[Y];
  1740.   Result.ImagPart[Z] := beta * QStart.ImagPart[Z] + t * QEnd.ImagPart[Z];
  1741.   Result.RealPart := beta * QStart.RealPart + t * QEnd.RealPart;
  1742. end;
  1743.  
  1744. //----------------------------------------------------------------------------------------------------------------------
  1745.  
  1746. // VectorAffineCombine
  1747. //
  1748. function VectorAffineCombine(const V1, V2: TAffineVector; F1, F2: Single): TAffineVector;
  1749. begin
  1750.   Result[X] := (F1 * V1[X]) + (F2 * V2[X]);
  1751.   Result[Y] := (F1 * V1[Y]) + (F2 * V2[Y]);
  1752.   Result[Z] := (F1 * V1[Z]) + (F2 * V2[Z]);
  1753. end;
  1754.  
  1755. // VectorAffineCombine3
  1756. //
  1757. function  VectorAffineCombine3(const V1, V2, V3: TAffineVector; F1, F2, F3: Single): TAffineVector;
  1758. begin
  1759.   Result[X] := (F1 * V1[X]) + (F2 * V2[X]) + (F3 * V3[X]);
  1760.   Result[Y] := (F1 * V1[Y]) + (F2 * V2[Y]) + (F3 * V3[Y]);
  1761.   Result[Z] := (F1 * V1[Z]) + (F2 * V2[Z]) + (F3 * V3[Z]);
  1762. end;
  1763.  
  1764. //----------------------------------------------------------------------------------------------------------------------
  1765.  
  1766. function VectorCombine(const V1, V2: TVector; F1, F2: Single): TVector;
  1767.  
  1768. // makes a linear combination of two vectors and return the result
  1769.  
  1770. begin
  1771.   Result[X] := (F1 * V1[X]) + (F2 * V2[X]);
  1772.   Result[Y] := (F1 * V1[Y]) + (F2 * V2[Y]);
  1773.   Result[Z] := (F1 * V1[Z]) + (F2 * V2[Z]);
  1774.   Result[W] := (F1 * V1[W]) + (F2 * V2[W]);
  1775. end;
  1776.  
  1777. //----------------------------------------------------------------------------------------------------------------------
  1778.  
  1779. function MatrixDecompose(const M: TMatrix; var Tran: TTransformations): Boolean; register;
  1780.  
  1781. // Author: Spencer W. Thomas, University of Michigan
  1782. //
  1783. // MatrixDecompose - Decompose a non-degenerated 4x4 transformation matrix into
  1784. // the sequence of transformations that produced it.
  1785. //
  1786. // The coefficient of each transformation is returned in the corresponding
  1787. // element of the vector Tran.
  1788. //
  1789. // Returns true upon success, false if the matrix is singular.
  1790.  
  1791. var I, J: Integer;
  1792.     LocMat,
  1793.     pmat,
  1794.     invpmat,
  1795.     tinvpmat: TMatrix;
  1796.     prhs,
  1797.     psol: TVector;
  1798.     Row: array[0..2] of TAffineVector;
  1799.  
  1800. begin
  1801.   Result := False;
  1802.   locmat := M;
  1803.   // normalize the matrix
  1804.   if locmat[W, W] = 0 then Exit;
  1805.   for I := 0 to 3 do
  1806.     for J := 0 to 3 do
  1807.       locmat[I, J] := locmat[I, J] / locmat[W, W];
  1808.  
  1809.   // pmat is used to solve for perspective, but it also provides
  1810.   // an easy way to test for singularity of the upper 3x3 component.
  1811.  
  1812.   pmat := locmat;
  1813.   for I := 0 to 2 do pmat[I, W] := 0;
  1814.   pmat[W, W] := 1;
  1815.  
  1816.   if MatrixDeterminant(pmat) = 0 then Exit;
  1817.  
  1818.   // First, isolate perspective.  This is the messiest.
  1819.   if (locmat[X, W] <> 0) or
  1820.      (locmat[Y, W] <> 0) or
  1821.      (locmat[Z, W] <> 0) then
  1822.   begin
  1823.     // prhs is the right hand side of the equation.
  1824.     prhs[X] := locmat[X, W];
  1825.     prhs[Y] := locmat[Y, W];
  1826.     prhs[Z] := locmat[Z, W];
  1827.     prhs[W] := locmat[W, W];
  1828.  
  1829.     // Solve the equation by inverting pmat and multiplying
  1830.     // prhs by the inverse.  (This is the easiest way, not
  1831.     // necessarily the best.)
  1832.  
  1833.     invpmat := pmat;
  1834.     MatrixInvert(invpmat);
  1835.     MatrixTranspose(invpmat);
  1836.     psol := VectorTransform(prhs, tinvpmat);
  1837.  
  1838.     // stuff the answer away
  1839.     Tran[ttPerspectiveX] := psol[X];
  1840.     Tran[ttPerspectiveY] := psol[Y];
  1841.     Tran[ttPerspectiveZ] := psol[Z];
  1842.     Tran[ttPerspectiveW] := psol[W];
  1843.  
  1844.     // clear the perspective partition
  1845.     locmat[X, W] := 0;
  1846.     locmat[Y, W] := 0;
  1847.     locmat[Z, W] := 0;
  1848.     locmat[W, W] := 1;
  1849.   end
  1850.   else
  1851.   begin
  1852.     // no perspective
  1853.     Tran[ttPerspectiveX] := 0;
  1854.     Tran[ttPerspectiveY] := 0;
  1855.     Tran[ttPerspectiveZ] := 0;
  1856.     Tran[ttPerspectiveW] := 0;
  1857.   end;
  1858.  
  1859.   // next take care of translation (easy)
  1860.   for I := 0 to 2 do
  1861.   begin
  1862.     Tran[TTransType(Ord(ttTranslateX) + I)] := locmat[W, I];
  1863.     locmat[W, I] := 0;
  1864.   end;
  1865.  
  1866.   // now get scale and shear
  1867.   for I := 0 to 2 do
  1868.   begin
  1869.     row[I, X] := locmat[I, X];
  1870.     row[I, Y] := locmat[I, Y];
  1871.     row[I, Z] := locmat[I, Z];
  1872.   end;
  1873.  
  1874.   // compute X scale factor and normalize first row
  1875.   Tran[ttScaleX] := Sqr(VectorNormalize(row[0])); // ml: calculation optimized
  1876.  
  1877.   // compute XY shear factor and make 2nd row orthogonal to 1st
  1878.   Tran[ttShearXY] := VectorAffineDotProduct(row[0], row[1]);
  1879.   row[1] := VectorAffineCombine(row[1], row[0], 1, -Tran[ttShearXY]);
  1880.  
  1881.   // now, compute Y scale and normalize 2nd row
  1882.   Tran[ttScaleY] := Sqr(VectorNormalize(row[1])); // ml: calculation optimized
  1883.   Tran[ttShearXY] := Tran[ttShearXY]/Tran[ttScaleY];
  1884.  
  1885.   // compute XZ and YZ shears, orthogonalize 3rd row
  1886.   Tran[ttShearXZ] := VectorAffineDotProduct(row[0], row[2]);
  1887.   row[2] := VectorAffineCombine(row[2], row[0], 1, -Tran[ttShearXZ]);
  1888.   Tran[ttShearYZ] := VectorAffineDotProduct(row[1], row[2]);
  1889.   row[2] := VectorAffineCombine(row[2], row[1], 1, -Tran[ttShearYZ]);
  1890.  
  1891.   // next, get Z scale and normalize 3rd row
  1892.   Tran[ttScaleZ] := Sqr(VectorNormalize(row[1])); // (ML) calc. optimized
  1893.   Tran[ttShearXZ] := Tran[ttShearXZ] / tran[ttScaleZ];
  1894.   Tran[ttShearYZ] := Tran[ttShearYZ] / Tran[ttScaleZ];
  1895.  
  1896.   // At this point, the matrix (in rows[]) is orthonormal.
  1897.   // Check for a coordinate system flip.  If the determinant
  1898.   // is -1, then negate the matrix and the scaling factors.
  1899.   if VectorAffineDotProduct(row[0], VectorCrossProduct(row[1], row[2])) < 0 then
  1900.     for I := 0 to 2 do
  1901.     begin
  1902.       Tran[TTransType(Ord(ttScaleX) + I)] := -Tran[TTransType(Ord(ttScaleX) + I)];
  1903.       row[I, X] := -row[I, X];
  1904.       row[I, Y] := -row[I, Y];
  1905.       row[I, Z] := -row[I, Z];
  1906.     end;
  1907.  
  1908.   // now, get the rotations out, as described in the gem
  1909.   Tran[ttRotateY] := arcsin(-row[0, Z]);
  1910.   if cos(Tran[ttRotateY]) <> 0 then
  1911.   begin
  1912.     Tran[ttRotateX] := arctan2(row[1, Z], row[2, Z]);
  1913.     Tran[ttRotateZ] := arctan2(row[0, Y], row[0, X]);
  1914.   end
  1915.   else
  1916.   begin
  1917.     tran[ttRotateX] := arctan2(row[1, X], row[1, Y]);
  1918.     tran[ttRotateZ] := 0;
  1919.   end;
  1920.   // All done!
  1921.   Result := True;
  1922. end;
  1923.  
  1924. //----------------------------------------------------------------------------------------------------------------------
  1925.  
  1926. function VectorDblToFlt(const V: THomogeneousDblVector): THomogeneousVector; assembler;
  1927.  
  1928. // converts a vector containing double sized values into a vector with single sized values
  1929.  
  1930. asm
  1931.               FLD  QWORD PTR [EAX]
  1932.               FSTP DWORD PTR [EDX]
  1933.               FLD  QWORD PTR [EAX + 8]
  1934.               FSTP DWORD PTR [EDX + 4]
  1935.               FLD  QWORD PTR [EAX + 16]
  1936.               FSTP DWORD PTR [EDX + 8]
  1937.               FLD  QWORD PTR [EAX + 24]
  1938.               FSTP DWORD PTR [EDX + 12]
  1939. end;
  1940.  
  1941. //----------------------------------------------------------------------------------------------------------------------
  1942.  
  1943. function VectorAffineDblToFlt(const V: TAffineDblVector): TAffineVector; assembler;
  1944.  
  1945. // converts a vector containing double sized values into a vector with single sized values
  1946.  
  1947. asm
  1948.               FLD  QWORD PTR [EAX]
  1949.               FSTP DWORD PTR [EDX]
  1950.               FLD  QWORD PTR [EAX + 8]
  1951.               FSTP DWORD PTR [EDX + 4]
  1952.               FLD  QWORD PTR [EAX + 16]
  1953.               FSTP DWORD PTR [EDX + 8]
  1954. end;
  1955.  
  1956. //----------------------------------------------------------------------------------------------------------------------
  1957.  
  1958. function VectorAffineFltToDbl(const V: TAffineVector): TAffineDblVector; assembler;
  1959.  
  1960. // converts a vector containing single sized values into a vector with double sized values
  1961.  
  1962. asm
  1963.               FLD  DWORD PTR [EAX]
  1964.               FSTP QWORD PTR [EDX]
  1965.               FLD  DWORD PTR [EAX + 8]
  1966.               FSTP QWORD PTR [EDX + 4]
  1967.               FLD  DWORD PTR [EAX + 16]
  1968.               FSTP QWORD PTR [EDX + 8]
  1969. end;
  1970.  
  1971. //----------------------------------------------------------------------------------------------------------------------
  1972.  
  1973. function VectorFltToDbl(const V: TVector): THomogeneousDblVector; assembler;
  1974.  
  1975. // converts a vector containing single sized values into a vector with double sized values
  1976.  
  1977. asm
  1978.               FLD  DWORD PTR [EAX]
  1979.               FSTP QWORD PTR [EDX]
  1980.               FLD  DWORD PTR [EAX + 8]
  1981.               FSTP QWORD PTR [EDX + 4]
  1982.               FLD  DWORD PTR [EAX + 16]
  1983.               FSTP QWORD PTR [EDX + 8]
  1984.               FLD  DWORD PTR [EAX + 24]
  1985.               FSTP QWORD PTR [EDX + 12]
  1986. end;
  1987.  
  1988. //----------------- coordinate system manipulation functions -----------------------------------------------------------
  1989.  
  1990. // Turn (Y axis)
  1991. //
  1992. function Turn(const Matrix: TMatrix; Angle: Single): TMatrix;
  1993. begin
  1994.    Result := MatrixMultiply(Matrix, CreateRotationMatrix(MakeAffineVector(Matrix[1]), Angle));
  1995. end;
  1996.  
  1997. // Turn (direction)
  1998. //
  1999. function Turn(const Matrix: TMatrix; const MasterUp: TAffineVector; Angle: Single): TMatrix;
  2000. begin
  2001.    Result := MatrixMultiply(Matrix, CreateRotationMatrix(MasterUp, Angle));
  2002. end;
  2003.  
  2004. // Pitch (X axis)
  2005. //
  2006. function Pitch(const Matrix: TMatrix; Angle: Single): TMatrix;
  2007. begin
  2008.    Result := MatrixMultiply(Matrix, CreateRotationMatrix(MakeAffineVector(Matrix[0]), Angle));
  2009. end;
  2010.  
  2011. // Pitch (direction)
  2012. //
  2013. function Pitch(const Matrix: TMatrix; const MasterRight: TAffineVector; Angle: Single): TMatrix; overload;
  2014. begin
  2015.    Result := MatrixMultiply(Matrix, CreateRotationMatrix(MasterRight, Angle));
  2016. end;
  2017.  
  2018. // Roll (Z axis)
  2019. //
  2020. function Roll(const Matrix: TMatrix; Angle: Single): TMatrix;
  2021. begin
  2022.    Result := MatrixMultiply(Matrix, CreateRotationMatrix(MakeAffineVector(Matrix[2]), Angle));
  2023. end;
  2024.  
  2025. // Roll (direction)
  2026. //
  2027. function Roll(const Matrix: TMatrix; const MasterDirection: TAffineVector; Angle: Single): TMatrix; overload;
  2028. begin
  2029.    Result := MatrixMultiply(Matrix, CreateRotationMatrix(MasterDirection, Angle));
  2030. end;
  2031.  
  2032. // MakeShadowMatrix
  2033. //
  2034. function MakeShadowMatrix(const planePoint, planeNormal, lightPos : TVector) : TMatrix;
  2035. var
  2036.    planeNormal3, dot : Single;
  2037. begin
  2038.     // Find the last coefficient by back substitutions
  2039.     planeNormal3:=-( planeNormal[0]*planePoint[0]
  2040.                    +planeNormal[1]*planePoint[1]
  2041.                    +planeNormal[2]*planePoint[2]);
  2042.     // Dot product of plane and light position
  2043.     dot:= planeNormal[0]*lightPos[0]
  2044.         +planeNormal[1]*lightPos[1]
  2045.         +planeNormal[2]*lightPos[2]
  2046.         +planeNormal3  *lightPos[3];
  2047.     // Now do the projection
  2048.     // First column
  2049.    Result[0][0] :=  dot - lightPos[0] * planeNormal[0];
  2050.    Result[1][0] :=      - lightPos[0] * planeNormal[1];
  2051.    Result[2][0] :=      - lightPos[0] * planeNormal[2];
  2052.    Result[3][0] :=      - lightPos[0] * planeNormal3;
  2053.     // Second column
  2054.     Result[0][1] :=      - lightPos[1] * planeNormal[0];
  2055.     Result[1][1] :=  dot - lightPos[1] * planeNormal[1];
  2056.     Result[2][1] :=      - lightPos[1] * planeNormal[2];
  2057.     Result[3][1] :=      - lightPos[1] * planeNormal3;
  2058.     // Third Column
  2059.     Result[0][2] :=      - lightPos[2] * planeNormal[0];
  2060.     Result[1][2] :=      - lightPos[2] * planeNormal[1];
  2061.     Result[2][2] :=  dot - lightPos[2] * planeNormal[2];
  2062.     Result[3][2] :=      - lightPos[2] * planeNormal3;
  2063.     // Fourth Column
  2064.     Result[0][3] :=      - lightPos[3] * planeNormal[0];
  2065.     Result[1][3] :=      - lightPos[3] * planeNormal[1];
  2066.     Result[2][3] :=      - lightPos[3] * planeNormal[2];
  2067.     Result[3][3] :=  dot - lightPos[3] * planeNormal3;
  2068. end;
  2069.  
  2070. end.
  2071.  
  2072.  
  2073.